//
//  ReadMapper.cpp
//  germ
//
//  Created by Vinicio Reynoso on 4/3/2010.
//  Purpose: class representing a kmer-decomposed sequence
//

#include "ReadMapper.h"
#include "TabFile.h"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <cstdlib>
#include <utility>
#include <stdint.h>

using namespace std;
//------------------------------------ constructor ---------------------------------------
ReadMapper::ReadMapper (int prefix_size):_prefix_size(prefix_size)
{
	_n_poss_prefixes = 0;
	_reset_data();
	
	// Calculate the number of possible prefixes.
	_n_poss_prefixes = 1;
	_n_poss_prefixes <<= _prefix_size;
	
	// Reserve space for 'binning' the sequence based on prefixes.
	_prsnt_prefixes = vBool(_n_poss_prefixes, false);
	_bins.reserve(_n_poss_prefixes);
	
	vector<unsigned int> * temp = NULL;
	for (u32 i = 0; i < _n_poss_prefixes; ++i) {
		_bins.push_back(temp);
	}
	
	// Make the prefix bit mask
	_prefix_bit_mask = 1;
	for (int i = 1; i < _prefix_size; ++i) {
		_prefix_bit_mask <<= 1;
		_prefix_bit_mask += 1;
	}
	
	// Reserve space for read mappings
	_mapping_ry.reserve(500000); _mapping_ws.reserve(500000);
	_mapping_mm.reserve(500000); _mapping_loc.reserve(500000);
	_mapping_gene.reserve(500000); _mapping_gene_strand.reserve(500000);
	_mapping_accession.reserve(500000); _mapping_org.reserve(500000);
}
//--------------------------------- destructor -------------------------------------------
ReadMapper::~ReadMapper()
{
	_reset_data();
	for (u32 i = 0; i < _n_poss_prefixes; ++i)
		delete _bins[i];
  
	_bins.clear();
}
//--------------------------------- _reset_data ------------------------------------------
void ReadMapper::_reset_data (void)
{
	_seqLoaded = false;
	_isAnnotated = false;
	_mappings_sorted = false;
	
	// Clear kmer information.
	ry_kmers.clear(); ws_kmers.clear(); locs.clear();
	
	// Clear protein information.
	protein_ranges.clear(); protein_syns.clear(); protein_names.clear();
	protein_eHits.clear(); protein_cHits.clear(); protein_products.clear();
	protein_strand.clear();
	
	// Clear RNA information.
	rna_ranges.clear(); rna_syns.clear(); rna_names.clear(); rna_eHits.clear();
	rna_cHits.clear(); rna_products.clear(); rna_strand.clear();
	
	// Clear read mapping information.
	_mapping_ry.clear(); _mapping_ws.clear(); _mapping_mm.clear(); _mapping_loc.clear();
	_mapping_gene.clear(); _mapping_gene_strand.clear(); _mapping_accession.clear();
	_mapping_org.clear(); _fmapping_ry.clear(); _fmapping_ws.clear();
	_fmapping_mm.clear(); _fmapping_loc.clear(); _fmapping_gene.clear();
	_fmapping_gene_strand.clear(); _fmapping_acession.clear(); _fmapping_org.clear();
  
	// Reset counts.
	_seq_length = _kmer_size = _n_kmers = _n_close = _n_fclose = _n_fexact =
	_n_exact_dup = _n_proteins = _n_exact_prot_hits = _n_close_prot_hits =
	_n_rna = _n_exact_rna_hits = _n_close_rna_hits = _n_eNCoding_hits =
	_n_cNCoding_hits = _t_spent_mapping = 0;
	
	// Clear bins.
	for (u32 i = 0; i < _n_poss_prefixes; ++i) {
		_prsnt_prefixes[i] = false;
    if (_bins[i] != NULL)
      _bins[i]->clear();
	}
	
	// Reset sequence information.
	_seq_source = "";
	_FASTA_file = "";
	_seq_files_Dir = "";
}
//--------------------------- reset_counts(void) -----------------------------------------
//should be used before mapping a new read-set to the sequence
void ReadMapper::reset_counts (void)
{
	// Reset gene counters.
	for (int i = 0; i < _n_proteins; ++i) {
		protein_eHits[i] = 0;
		protein_cHits[i] = 0;
	}
	for (int i = 0; i < _n_rna; ++i) {
		rna_eHits[i] = 0;
		rna_cHits[i] = 0;
	}
	
	// Clear mappings
	_mapping_ry.clear(); _mapping_ws.clear(); _mapping_mm.clear(); _mapping_loc.clear();
	_mapping_gene.clear(); _mapping_gene_strand.clear(); _mapping_accession.clear();
	_mapping_org.clear(); _fmapping_ry.clear(); _fmapping_ws.clear();
	_fmapping_mm.clear(); _fmapping_loc.clear(); _fmapping_gene.clear();
	_fmapping_gene_strand.clear(); _fmapping_acession.clear(); _fmapping_org.clear();
	
	// Reset counts
	_n_close = _n_fclose = _n_fexact = _n_exact_dup = _n_exact_prot_hits =
	_n_close_prot_hits = _n_exact_rna_hits = _n_close_rna_hits =
	_n_eNCoding_hits = _n_cNCoding_hits = _t_spent_mapping = 0;
	_mappings_sorted = false;
}
//---------------------------- load_seq() ------------------------------------------------
void ReadMapper::load_seq (int kmersize,
                           string fna_file,
                           string seqFileDir,
                           bool annotate)
{
	if (kmersize > 64) {
		cerr << "IN: ReadMapper::load_seq(...);" << endl;
		cerr << "\tERROR: Cannot represent a kmer longer than 64 nts. _kmer_size= ";
		cerr << _kmer_size << endl;
		_seqLoaded = false;
		return;
	}
	if (kmersize < _prefix_size) {
		cerr << "IN: ReadMapper::load_seq(...);" << endl;
		cerr << "\tERROR: kmer size is smaller than the ry-prefix size" << endl;
		_seqLoaded = false;
		return;
	}
	_reset_data();
	_kmer_size = kmersize;
	_FASTA_file = fna_file;
	_seq_files_Dir = seqFileDir;
	
	// Get the accession number.
	_accession = _FASTA_file.substr(0, _FASTA_file.find_first_of("."));
	
	string gbk_file;
	if(_seq_files_Dir == ""){
		gbk_file = _accession + ".gbk";
	}else{
		gbk_file = _seq_files_Dir + _accession + ".gbk";
		_FASTA_file = _seq_files_Dir + _FASTA_file;
	}
  
	
	/*  ---------------------------- gbk file example ----------------------------------
	 *  LOCUS       NS_000196             243619 bp    DNA     linear   CON 13-APR-2009
	 *	DEFINITION  Borrelia garinii pBi plasmid pseudoscaffold.
	 *	ACCESSION   NS_000196
	 *	VERSION     NS_000196.1  GI:226526126
	 *	KEYWORDS    .
	 */
	
	// Determine if the genome is circular.
	TabFile gbkTF;
	if (!gbkTF.open(gbk_file, false)) {
		cerr << "IN: ReadMapper::load_seq(...);" << endl;
		cerr << "\tERROR: could not open file= " << gbk_file << endl;
	} else {
		vector<string> buff;
		gbkTF.getField(buff);
		
		string line = buff[0];
		size_t found_cir = line.find("circular");
		
		if(found_cir != string::npos)
			_isCircular = true;
		else
			_isCircular = false;
		
		gbkTF.close();
	}
	
	// Load the sequence.
	vector<char> seq;
	
	TabFile fnaTF;
	if (!fnaTF.open(_FASTA_file)) {
		cerr << "IN: ReadMapper::load_seq(...);" << endl;
		cerr << "\tERROR: could not open file= " << _FASTA_file << endl;
	} else {
    
		int n_rows = fnaTF.nRows(); fnaTF.fp_data_start();
		seq.reserve(n_rows*70);
		
		_fasta_header = fnaTF.getHeaderString();
		_seq_source = _fasta_header.substr(_fasta_header.find_last_of("|")+2);
    
		vector<string> buff;
		string::iterator i;
		
		while (fnaTF.getField(buff)) {
			if (buff.empty()) continue;
			
			for (i = buff[0].begin(); i != buff[0].end(); ++i)
				seq.push_back(*i);
		}
		fnaTF.close();
	}
	
	// Get the number of k-mers.
	_seq_length = seq.size();
	int nKmers;
	if(_isCircular)
		nKmers = _seq_length * 2;
	else
		nKmers = (_seq_length - _kmer_size + 1) * 2;
  
	// Reserve space for the ry, ws, and their locs
	ry_kmers.reserve(nKmers);
	ws_kmers.reserve(nKmers);
	locs.reserve(nKmers);
	
	// If the sequence is circular, add (_kmer_size - 1) beginning nucleotides to the end
	int last_kmer_index;
	if (_isCircular) {
		for (int i = 0; i < _kmer_size-1; ++i)
			seq.push_back(seq[i]);
    
		last_kmer_index = _seq_length - 1;
	} else {
		last_kmer_index = _seq_length - _kmer_size;
	}
	
	// Convert the sequence into ry-ws integer pairs, which represent a single k-mer
	bool found_bad_char;
	u64 ry_int, ws_int, rc_ry_int, rc_ws_int, on;
	
	for (int seq_index = 0; seq_index <= last_kmer_index; ++seq_index) {
		ry_int = ws_int = rc_ry_int = rc_ws_int = 0;
		on = 1;
		
		int i = seq_index;
		found_bad_char = false;
		
		for (int j = 0; j < _kmer_size; ++j) {	// Kmer loop
			switch (seq[i]) {
				case 'A' :
					// Add "A" (00)
					ry_int <<= 1;	ws_int <<= 1;
					
					// Add "T" (10) to reverse complement
					rc_ry_int += on;
					
					on <<= 1;
					break;
					
				case 'G' :
					// Add "G" (01)
					ry_int <<= 1;
					ws_int <<= 1;	ws_int += 1;
					
					// Add "C" (11) to reverse complement
					rc_ry_int += on;
					rc_ws_int += on;
					
					on <<= 1;
					break;
					
				case 'T' :
					// Add "T" (10)
					ry_int <<= 1;	ry_int += 1;
					ws_int <<= 1;
					
					on <<= 1;
					break;
					
				case 'C' :
					//Add "C" (11) to ints
					ry_int <<= 1;	ry_int += 1;
					ws_int <<= 1;	ws_int += 1;
					
					//Add "G" (01) to rc ints
					rc_ws_int += on;
					
					on <<= 1;
					break;
					
				default :
					found_bad_char = true;
			}
			if(found_bad_char) break;
			++i;
		}
		
		if (found_bad_char) {
			seq_index = i;
			continue;
		}
		// Record k-mers
		ry_kmers.push_back(ry_int); ry_kmers.push_back(rc_ry_int);
		ws_kmers.push_back(ws_int); ws_kmers.push_back(rc_ws_int);
		
		locs.push_back(seq_index + 1);		// Sequence doesn't start at zero so add one.
		locs.push_back((seq_index + 1) * -1);	// Negative loc designates reverse
    // complement.
	}
	seq.clear();
	
	// Sort ry-ws integer pairs based on ry-integers.
	_n_kmers = ry_kmers.size();
	_sort_kmer_int_pairs(0, _n_kmers-1, true, false);
	
	// Now sort based on ws-integers
	int left, right;
	int end = _n_kmers - 1;
	for (left = 0; left < end; ++left) {
		right = left+1;
		
		if (ry_kmers[left] == ry_kmers[right]) {
      
			while ( (right+1 < _n_kmers) && (ry_kmers[left] == ry_kmers[right+1]) )
				++right;
      
			_sort_kmer_int_pairs(left, right, false, true);
			
			left = right;
		}
	}
	// Annotate the sequence
	if(annotate)
		_annotate_seq();
  
	// Fill bins with 1-transversion redundancy for mapping close reads
	u32 prefix, prefix_w_tv, bit_twiddle;
	
	for (int i = 0; i < _n_kmers; ++i) {
		prefix = _prefix_bit_mask & ry_kmers[i];
		
		_prsnt_prefixes[prefix] = true;
    
    if (_bins[prefix] == NULL)
      _bins[prefix] = new vector<unsigned>();
    
		_bins[prefix]->push_back((unsigned)i);
	}
	
	// Get the size of each bin.
	_bin_sizes = vector<int>(_n_poss_prefixes, 0);
	for(u32 i = 0; i < _n_poss_prefixes; ++i) {
    if (_bins[i] == NULL)
      _bin_sizes[i] = 0;
    else
      _bin_sizes[i] = _bins[i]->size();
  }
  
	// Add 1-transversion redundancy to the bins
	for (u32 i = 0; i < _n_poss_prefixes; ++i) {
		if (_prsnt_prefixes[i] == true) {
			
			for (int j = 0; j < _prefix_size; ++j) {
				bit_twiddle = 1;
				bit_twiddle <<= j;
				
				prefix_w_tv = bit_twiddle ^ i;
				
				if (_prsnt_prefixes[prefix_w_tv] == true) {
					
					int n_to_copy = _bin_sizes[prefix_w_tv];
					for (int k = 0; k < n_to_copy; ++k)
						_bins[i]->push_back(_bins[prefix_w_tv]->at(k));
				}
			}
		}
	}
	_seqLoaded = true;
}
//----------------------------- _sort_kmer_int_pairs() -----------------------------------
void ReadMapper::_sort_kmer_int_pairs (int left, int right, bool ry, bool ws)
{
	if ( ry && ws ) {
		cerr << "IN: ReadMapper::_sort_kmer_int_pairs(...);" << endl;
		cerr << "\tERROR: Cannot sort as such\t ry: " << ry << " ws: " << ws << endl;
		return;
	}
	
  int i = left;
  int j = right;
  int temp_loc;
  
  u64 temp;
  u64 x = 0;
  
  if(ry) x = ry_kmers[(left+right)/2];
  if(ws) x = ws_kmers[(left+right)/2];
  
  do{
    if (ry) {
      
      while ( (ry_kmers[i] < x) && (i < right) )    i++;
      
      while ( (x < ry_kmers[j]) && (j > left) )     j--;
    	
    } else {
      
      while ( (ws_kmers[i] < x) && (i < right) )    i++;
      
      while ( (x < ws_kmers[j]) && (j > left) )     j--;
    }
    if (i <= j) {
      temp = ry_kmers[i];	// Swap ry integers
      ry_kmers[i] = ry_kmers[j];
      ry_kmers[j] = temp;
      
      temp = ws_kmers[i];	// Swap ws integers at the same indices
      ws_kmers[i] = ws_kmers[j];
      ws_kmers[j] = temp;
      
      temp_loc = locs[i];	// Swap locs at the same indices
      locs[i] = locs[j];
      locs[j] = temp_loc;
      i++; j--;
    }
    
  } while (i <= j);
  
  if(i < right)
    _sort_kmer_int_pairs(i, right, ry, ws);
  
  if(left < j)
    _sort_kmer_int_pairs(left, j, ry, ws);
}
//----------------------------------- _annotate_seq() ------------------------------------
void ReadMapper::_annotate_seq (void)
{
	if (_isAnnotated) {
		cerr << "IN: ReadMapper::_annotate_seq(void);" << endl;
		cerr << "\tERROR: sequence has already been annotated" << _seq_source << endl;
		return;
	}
	
	if (_seq_files_Dir == "") {
		_ptt_file = _accession + ".ptt";
		_rnt_file = _accession + ".rnt";
	} else {
		_ptt_file = _seq_files_Dir + _accession + ".ptt";
		_rnt_file = _seq_files_Dir + _accession + ".rnt";
	}
	
	vector<string> buff;
	
	/*  ---------------------------- .ptt example ---------------------------------
	 *  Borrelia garinii pBi plasmid pseudoscaffold - 1..243619
	 *  338 proteins
	 *  Location	Strand	Length	PID	Gene	Synonym	Code	COG	Product
	 *  427..1122	-	231	56560866	-	BGP001	-	-	hypothetical protein
	 *  1181..1732	-	183	56560867	-	BGP002  -	-	hypothetical protein
	 */
  
	TabFile pttTF;
	if (!pttTF.open(_ptt_file, false)) {
		// Assuming all sequences have a .ptt file
	} else {
		// Get the number of proteins
		_n_proteins = pttTF.nRows() - 3;
		pttTF.fp_data_start();
		
		// Reserve space for protein information
		protein_ranges.reserve(_n_proteins);
		protein_names.reserve(_n_proteins);
		protein_syns.reserve(_n_proteins);
		protein_products.reserve(_n_proteins);
		protein_strand.reserve(_n_proteins);
		
		protein_eHits = vector<int>(_n_proteins, 0);
		protein_cHits = vector<int>(_n_proteins, 0);
		
		// Get the protein information
		pttTF.getNextRow(buff, 3);	//advance file pointer to beginning of protein info
		
		vector<int> range;
		while (pttTF.getField(buff)) {
			if (buff.empty()) continue;
			
			// Get the gene range (start..end)
			string gene_range = buff[0];
			string start = gene_range.substr(0, gene_range.find_first_of("."));
			string end = gene_range.substr(gene_range.find_first_of(".")+2);
			range.push_back(atoi(start.c_str()));
			range.push_back(atoi(end.c_str()));
			protein_ranges.push_back(range);
			range.clear();
			
			// Get the coding strand of the gene
			char strand = buff[1].at(0);
			protein_strand.push_back(strand);
			protein_names.push_back(buff[4]);	// Get the gene name
			protein_syns.push_back(buff[5]);	// Get the synonym
			protein_products.push_back(buff[8]);	// Get the protein product
		}
		pttTF.close();
	}
	
	/*  ---------------------------- .rnt example -------------------------------------
	 *  Uncultured Termite group 1 bacterium phylotype Rs-D17, complete genome - 1..112
	 *	48 RNAs
	 *	Location	Strand	Length	PID	Gene	Synonym	Code	COG	Product
	 *	23299..23374	+	76	189485003	-	TGRD_tRNA01	-	-	Val tRNA
	 *	34638..34727    +	90	189485003	-	TGRD_tRNA02	-	-	Ser tRNA
	 */
	
	TabFile rntTF;
	if (!rntTF.open(_rnt_file, false)) {
		// Assuming all sequences have a .rnt file
	} else {
		// Get the number of RNA genes
		_n_rna = rntTF.nRows() - 3;
		rntTF.fp_data_start();
		
		// Reserve space
		rna_ranges.reserve(_n_rna);
		rna_eHits = vector<int>(_n_rna,0);
		rna_cHits = vector<int>(_n_rna,0);
		rna_names.reserve(_n_rna);
		rna_syns.reserve(_n_rna);
		rna_products.reserve(_n_rna);
		rna_strand.reserve(_n_rna);
		
		// Get the RNA information
		rntTF.getNextRow(buff, 3);	// Advancing file pointer to beginning of rna info
    
		vector<int> range;
		while (rntTF.getField(buff)) {
			if (buff.empty()) continue;
      
			string gene_range = buff[0];	// Get the gene range
			string start = gene_range.substr(0, gene_range.find_first_of("."));
			string end = gene_range.substr(gene_range.find_first_of(".")+2);
			range.push_back(atoi(start.c_str()));
			range.push_back(atoi(end.c_str()));
			rna_ranges.push_back(range);
			range.clear();
			
			char strand = buff[1].at(0);	// Get the coding strand of the gene
			rna_strand.push_back(strand);
			rna_names.push_back(buff[4]);	// Get the gene name
			rna_syns.push_back(buff[5]);	// Get the synonym
			rna_products.push_back(buff[8]);	// Get the rna product
		}
		rntTF.close();
	}
	_isAnnotated = true;
}
//----------------------------------- map_reads() ----------------------------------------
// It's assumed that read ry-ws integer pairs are sorted and all duplicates removed.
void ReadMapper::map_reads (v64& ry_reads,
                            v64& ws_reads,
                            vBool& mapped)
{
	if (!_seqLoaded) return;
	
	int k = 0;  // Sequence k-mer index
	long int r = 0;  // Read index
	
	long int n_reads = ry_reads.size();
	vector<string> genes_hit;
	time_t start_t, end_t;
	
	// Map the reads
	start_t = time(NULL);
	while ( (r < n_reads) && (k < _n_kmers) ) {
		if (ry_reads[r] == ry_kmers[k]) {
			if (ws_reads[r] == ws_kmers[k]) {
        
				// The two ry-ws integer pairs are identical
				_record_mapping(ry_reads[r], ws_reads[r], locs[k], 0);
				mapped[r] = true;
				++_n_fexact;
				
				// Look for duplicate k-mers in the sequence
				while ( (k+1 < _n_kmers) && (ry_kmers[k] == ry_kmers[k+1]) ) {
					
					if (ws_kmers[k] == ws_kmers[k+1]) {
						++k;
						_record_mapping(ry_reads[r], ws_reads[r], locs[k], 0);
						++_n_exact_dup;
					} else {
						break;
					}
				}
				++r; ++k;	// Advance to the next read and sequence k-mer
			} else {
				// ry's were the same but not ws's.
				// Increment the index of the smaller ws int
				if(ws_reads[r] < ws_kmers[k])
					++r;
				else
					++k;
			}
		} else {
			// ry's were not the same. Increment the index of the smaller ry int
			if(ry_reads[r] < ry_kmers[k])
				++r;
			else
				++k;
		}
	}
	end_t = time(NULL);
	_t_spent_mapping += difftime(end_t, start_t);
}
//----------------------------------- map_close_reads() ----------------------------------
void ReadMapper::map_close_reads (v64& ry_reads,
                                  v64& ws_reads,
                                  vBool& mapped,
                                  int threshold)
{
	if (!_seqLoaded) return;
	
	// Map close reads
	u32 prefix;
	long int n_reads = ry_reads.size();
	
	time_t start_t = time(NULL);
	for (long int r = 0; r < n_reads; ++r) {
		if (mapped[r] == false) {
			prefix = _prefix_bit_mask & ry_reads[r];
			
			if (_check_bins(ry_reads[r], ws_reads[r], prefix, threshold))
				mapped[r] = true;
		}
	}
	
	time_t end_t = time(NULL);
	_t_spent_mapping += difftime(end_t, start_t);
}
//------------------------------------ _check_bins() -------------------------------------
// Checks all bins that are one transversion away from the prefix of the read's
// ry integer, in addition to the exact prefix, for close matches.
bool ReadMapper::_check_bins (u64 ry_read,
                              u64 ws_read,
                              u32 prefix,
                              int threshold)
{
	u64 x, y, z;
	u32 prefix_w_tv, bit_twiddle;
	bool found_close = false;
  vector<unsigned int>::iterator index_itr, end_of_bin;
  
  bit_twiddle = 0;
  int nBins = _prefix_size + 1;
  int mm;	// Number of mismatches.
  
  for (int t = 0; t < nBins; ++t) {
    prefix_w_tv = bit_twiddle ^ prefix;
    
    if (_prsnt_prefixes[prefix_w_tv] == true) {
      index_itr = _bins[prefix_w_tv]->begin();
			end_of_bin = _bins[prefix_w_tv]->end();
			
			do{
        
				x = ry_read ^ ry_kmers[*index_itr];
				
				mm = __builtin_popcountl(x);
				if (mm <= threshold) {
          
					y = x ^ (ws_read ^ ws_kmers[*index_itr]);
					z = y & ~x;
					
					mm += __builtin_popcountl(z);
					if (mm <= threshold) {
						_record_mapping(ry_read, ws_read, locs[*index_itr], mm);
						++_n_close;
						found_close = true;
					}
				}
				++index_itr;
			}while (index_itr != end_of_bin);
		}
		
		// Set next transversion position
		bit_twiddle = 1;
		bit_twiddle <<= t;
  }
	return found_close;
}
//--------------------------------- _record_mapping() ------------------------------------
//header: "ry_read\tid_read\tmm\tloc\tgene_syn\tgene_strand\taccession\torg\n"
void ReadMapper::_record_mapping (u64 ry_read,
                                  u64 ws_read,
                                  int loc,
                                  int mm)
{
	// Is it a close mapping?
	bool close;
	if(mm == 0)
		close = false;
	else
		close = true;
	
	vector<string> genes_hit;
	vector<char> strand_of_gene_hit;
	_genes_hit(loc, genes_hit, strand_of_gene_hit, close);
	
	int n_genes_hit = genes_hit.size();
	for (int i = 0; i < n_genes_hit; ++i) {
		_mapping_ry.push_back(ry_read);
		_mapping_ws.push_back(ws_read);
		_mapping_mm.push_back(mm);
		_mapping_loc.push_back(loc);
		_mapping_gene.push_back(genes_hit[i]);
		_mapping_gene_strand.push_back(strand_of_gene_hit[i]);
		_mapping_accession.push_back(_accession);
		_mapping_org.push_back(_seq_source);
	}
}
//------------------------------------- _genes_hit() -------------------------------------
// Gene hit counters are incremented only when the read maps onto the coding strand of
// the gene. However, regardless of which strand of the gene the read maps onto, the
// read's mapping information will still say that it mapped onto that particular gene.
void ReadMapper::_genes_hit (int loc,
                             vector<string>& genes_hit,
                             vector<char>& strand_of_gene_hit,
                             bool close)
{
	if(!_isAnnotated) return;
	genes_hit.clear();
	strand_of_gene_hit.clear();
	
	// Check whether the read mapped onto the forward (+) or reverse (-) strand
	char strand;
	if(loc < 0){
		strand = '-';
		loc *= -1;
	} else {
		strand = '+';
	}
	
	// Check if the read mapped onto a protein gene
	int start, end;
	for (int i = 0; i < _n_proteins; ++i) {
		start = protein_ranges[i][0] - (_kmer_size-1);
		end = protein_ranges[i][1] - 1;
		
		if (loc >= start) {
			if (loc <= end) {
				if (protein_strand[i] == strand) {
					// The read mapped onto the coding strand of the gene,
					// increment hit counter
					if (close)
						++protein_cHits[i];
					else
						++protein_eHits[i];
          
					genes_hit.push_back(protein_syns[i]);
					strand_of_gene_hit.push_back(protein_strand[i]);
				} else {
					// Read mapped onto the non-coding strand, so do not increment
					// hit counters
					genes_hit.push_back(protein_syns[i]);
					strand_of_gene_hit.push_back(protein_strand[i]);
				}
			}
		} else break;	// Since genes are sorted based on range, if the mapping loc is
    // smaller than the start (minus the kmer size) of any gene,
    // then the read did not map onto any of the other genes.
	}
	
	// Check if the read mapped onto a rna gene
	for (int i = 0; i < _n_rna; ++i) {
		start = rna_ranges[i][0] - (_kmer_size-1);
		end = rna_ranges[i][1] - 1;
		
		if (loc >= start) {
			if (loc <= end) {
				if (rna_strand[i] == strand) {
					// The read mapped onto the coding region of the rna gene
					if (close)
						++rna_cHits[i];
					else
						++rna_eHits[i];
          
					genes_hit.push_back(rna_syns[i]);
					strand_of_gene_hit.push_back(rna_strand[i]);
				} else {
					// Read mapped onto the non-coding strand, so do not increment
					// hit counters
					genes_hit.push_back(rna_syns[i]);
					strand_of_gene_hit.push_back(rna_strand[i]);
				}
			}
		}else break;
	}
	
	if (genes_hit.empty()) {	// Read mapped onto a non-coding region
		genes_hit.push_back("-");
		strand_of_gene_hit.push_back('*');
	}
}
//--------------------------- filter_mappings() - overloaded -----------------------------
void ReadMapper::filter_mappings (void)
{
	// Reserve space for filtered mappings
	long int n_mappings = _mapping_ry.size();
	_fmapping_ry.reserve(n_mappings); _fmapping_ws.reserve(n_mappings);
	_fmapping_mm.reserve(n_mappings); _fmapping_loc.reserve(n_mappings);
	_fmapping_gene.reserve(n_mappings);	_fmapping_gene_strand.reserve(n_mappings);
	_fmapping_acession.reserve(n_mappings); _fmapping_org.reserve(n_mappings);
	--n_mappings;
	
	// Sort the mappings
	if (!_mappings_sorted)
		sort_mappings();
	
	// Write the closest mapping for each read in each organism.
	long int left, right;
	for (left = 0; left < n_mappings; ++left) {
		right = left+1;
		
		if ( (_mapping_ry[left] == _mapping_ry[right]) &&
        (_mapping_ws[left] == _mapping_ws[right]) ) {
			
			while ( (right+1 < n_mappings) &&
             (_mapping_ry[left] == _mapping_ry[right+1]) &&
             (_mapping_ws[left] == _mapping_ws[right+1]) )
				++right;
      
			_record_best_mappings(left, right);
			left = right;
		} else {
			// The read only mapped once, record its mapping info
			_record_best_mappings(left, left);
		}
	}
}

// This version loads mappings from file, filters them, and then prints them back to file
// CAUTION: Any present mappings are erased!
void ReadMapper::filter_mappings (string& file)
{
	// Load the read mappings from file
	load_mappings(file, false);		// Assuming mappings are not sorted
  
	// Clear the read mappings from file
	TabFile read_mappingsTF;
	if (!read_mappingsTF.open(file)) {
		cerr << "IN: ReadMapper::filter_mappings(string& file);" << endl;
		cerr << "\tERROR: could not open file= " << file << endl;
		return;
	} else {
		vector<string> buff;
		buff = read_mappingsTF.getHeader();
		read_mappingsTF.clearDBFile();
		read_mappingsTF.setHeader(buff);
		read_mappingsTF.writeHeader();
		read_mappingsTF.close();
	}
  
	// Filter all the read mappings
	filter_mappings();
	write_mappings(file, true);
}
//------------------ _record_best_mappings(long int left, long int right) ----------------
// Assumes that at indices [left,right] are all the mappings for one particular read
// and that the mappings are sorted based on mismatch number
void ReadMapper::_record_best_mappings (long int left, long int right)
{
	// Record the first mapping since it's the closest or one of the closest
	_fmapping_ry.push_back(_mapping_ry[left]);
	_fmapping_ws.push_back(_mapping_ws[left]);
	_fmapping_mm.push_back(_mapping_mm[left]);
	
	if (_mapping_mm[left] > 0)
		++_n_fclose;
	
	_fmapping_loc.push_back(_mapping_loc[left]);
	_fmapping_gene.push_back(_mapping_gene[left]);
	_fmapping_gene_strand.push_back(_mapping_gene_strand[left]);
	_fmapping_acession.push_back(_mapping_accession[left]);
	_fmapping_org.push_back(_mapping_org[left]);
  
	if (left == right)	// Only one read mapping to record
		return;
	
	set<string> orgs_that_read_mapped_to;
	pair< set<string>::iterator, bool > insrt_res;
	orgs_that_read_mapped_to.insert(_mapping_accession[left]);
	int smallest_mm_n = _mapping_mm[left];
	
	for (long int i = left+1; i <= right; ++i) {
		// Note: mappings at indices left to right have been sorted based on
		// mismatch number
		if (_mapping_mm[i] == smallest_mm_n) {
			insrt_res = orgs_that_read_mapped_to.insert(_mapping_accession[i]);
			
			if (insrt_res.second == true) {
				// Read mapped just as closely in a different org.
				// Record that mapping as well.
				_fmapping_ry.push_back(_mapping_ry[i]);
				_fmapping_ws.push_back(_mapping_ws[i]);
				_fmapping_mm.push_back(_mapping_mm[i]);
				
				if (_mapping_mm[i] > 0)
					++_n_fclose;
				
				_fmapping_loc.push_back(_mapping_loc[i]);
				_fmapping_gene.push_back(_mapping_gene[i]);
				_fmapping_gene_strand.push_back(_mapping_gene_strand[i]);
				_fmapping_acession.push_back(_mapping_accession[i]);
				_fmapping_org.push_back(_mapping_org[i]);
			}
		}else return;	// The rest of the mappings were at a further distance.
    // (these get ignored)
	}
}
//----------------------------------- load_mappings(...) ---------------------------------
// This method will clear all mapping information before it loads from file
void ReadMapper::load_mappings (string& file, bool sorted)
{
	// Clear mapping information.
	_mapping_ry.clear(); _mapping_ws.clear(); _mapping_mm.clear(); _mapping_loc.clear();
	_mapping_gene.clear(); _mapping_gene_strand.clear(); _mapping_accession.clear();
	_mapping_org.clear(); _fmapping_ry.clear(); _fmapping_ws.clear();
	_fmapping_mm.clear(); _fmapping_loc.clear(); _fmapping_gene.clear();
	_fmapping_gene_strand.clear(); _fmapping_acession.clear(); _fmapping_org.clear();
	
	// Open the mapping file
	TabFile read_mappingsTF;
	if (!read_mappingsTF.open(file)) {
		cerr << "IN: ReadMapper::filter_mappings(string& file);" << endl;
		cerr << "\tERROR: could not open file= " << file << endl;
		return;
	} else {
		// Reserve space for the unfiltered mappings
		u64 n_mappings = read_mappingsTF.nRows();
		read_mappingsTF.fp_data_start();
		
		_mapping_ry.reserve(n_mappings); _mapping_ws.reserve(n_mappings);
		_mapping_mm.reserve(n_mappings); _mapping_loc.reserve(n_mappings);
		_mapping_gene.reserve(n_mappings); _mapping_gene_strand.reserve(n_mappings);
		_mapping_accession.reserve(n_mappings); _mapping_org.reserve(n_mappings);
		
		// Load all the mapping information.
		char strand;
		vector<string> buff;
		while (read_mappingsTF.getField(buff)) {
			if (buff.empty()) continue;
			
			_mapping_ry.push_back( strtoul(buff[0].c_str(),NULL,10) );
			_mapping_ws.push_back( strtoul(buff[1].c_str(),NULL,10) );
			_mapping_mm.push_back( atoi(buff[2].c_str()) );
			_mapping_loc.push_back( atoi(buff[3].c_str()) );
			_mapping_gene.push_back(buff[4]);
			strand = buff[5].at(0);
			_mapping_gene_strand.push_back(strand);
			_mapping_accession.push_back(buff[6]);
			_mapping_org.push_back(buff[7]);
		}
		read_mappingsTF.close();
	}
	if (sorted)
		_mappings_sorted = true;
	else
		_mappings_sorted = false;
}
//------------------------------- sort_mappings(void) ------------------------------------
void ReadMapper::sort_mappings (void)
{
	if (_mappings_sorted) return;
	
	// Sort mappings based on ry
	long int n_mappings = _mapping_ry.size();
	_sort_mappings(0, n_mappings-1, true, false, false);
	
	// Sort mappings based on ws
	long int left, right;
	for (left = 0; left < n_mappings-1; ++left) {
		right = left+1;
		
		if (_mapping_ry[left] == _mapping_ry[right]) {
			
			while ( (right+1 < n_mappings) &&
             (_mapping_ry[left] == _mapping_ry[right+1]) )
				++right;
      
			_sort_mappings(left, right, false, true, false);
			left = right;
		}
	}
	
	// Sort mappings based on mm number. Here, all the mappings for one particular
	// read get sorted.
	for (left = 0; left < n_mappings-1; ++left) {
		right = left+1;
		
		if ( (_mapping_ry[left] == _mapping_ry[right]) &&
        (_mapping_ws[left] == _mapping_ws[right]) ) {
      
			while ( (right+1 < n_mappings) &&
             (_mapping_ry[left] == _mapping_ry[right+1]) &&
             (_mapping_ws[left] == _mapping_ws[right+1]) )
				++right;
      
			_sort_mappings(left, right, false, false, true);
			left = right;
		}
	}
	_mappings_sorted = true;
}
//------------------------------ _sort_mappings(...) -------------------------------------
void ReadMapper::_sort_mappings (long int left,
                                 long int right,
                                 bool ry,
                                 bool ws,
                                 bool mm)
{
	// Can only sort on one thing at a time
	if ( (ry||ws) && mm ) {
		cerr << "IN: ReadMapper::_sort_mappings(...);" << endl;
		cerr << "\tERROR: Cannot sort as such\t ry: " << ry << " ws: " << ws;
		cerr << " mm: " << mm << endl;
		return;
	}
	
	long int i = left;
	long int j = right;
	
	int temp_int;
	int xx = 0;
  u64 temp_ul;
  u64 x = 0;
  
  string temp_str;
  char temp_char;
  
  if (ry)
    x = _mapping_ry[(left+right)/2];
  if (ws)
    x = _mapping_ws[(left+right)/2];
  if (mm)
    xx = _mapping_mm[(left+right)/2];
  
  do{
    if (ry) {
      while ( (_mapping_ry[i] < x) && (i < right) )    i++;
      while ( (x < _mapping_ry[j]) && (j > left) )     j--;
    }
    
    if (ws) {
      while ( (_mapping_ws[i] < x) && (i < right) )    i++;
      while ( (x < _mapping_ws[j]) && (j > left) )     j--;
    }
    
    if (mm) {
      while ( (_mapping_mm[i] < xx) && (i < right) )    i++;
      while ( (xx < _mapping_mm[j]) && (j > left) )     j--;
    }
    
    if (i <= j) {
      temp_ul = _mapping_ry[i];	// Swap ry integers
      _mapping_ry[i] = _mapping_ry[j];
      _mapping_ry[j] = temp_ul;
      
      temp_ul = _mapping_ws[i];	// Swap ws integers
      _mapping_ws[i] = _mapping_ws[j];
      _mapping_ws[j] = temp_ul;
      
      temp_int = _mapping_mm[i];	// Swap mismatch numbers
      _mapping_mm[i] = _mapping_mm[j];
      _mapping_mm[j] = temp_int;
      
      temp_int = _mapping_loc[i];	// Swap locs
      _mapping_loc[i] = _mapping_loc[j];
      _mapping_loc[j] = temp_int;
      
      temp_str = _mapping_gene[i];	// Swap gene synonyms
      _mapping_gene[i] = _mapping_gene[j];
      _mapping_gene[j] = temp_str;
      
			temp_char = _mapping_gene_strand[i];	// Swap gene strands
			_mapping_gene_strand[i] = _mapping_gene_strand[j];
			_mapping_gene_strand[j] = temp_char;
			
      temp_str = _mapping_accession[i];	// Swap accessions
      _mapping_accession[i] = _mapping_accession[j];
      _mapping_accession[j] = temp_str;
      
      temp_str = _mapping_org[i];			// Swap org names
      _mapping_org[i] = _mapping_org[j];
      _mapping_org[j] = temp_str;
      
      i++; j--;
    }
  } while (i <= j);
  
  if(i < right)
    _sort_mappings(i, right, ry, ws, mm);
  if(left < j)
    _sort_mappings(left, j, ry, ws, mm);
}
//----------------------------- write_mappings() -----------------------------------------
void ReadMapper::write_mappings (string& file, bool filtered)
{
	ofstream out(file.c_str(), ios_base::app);
	if (!out) {
		cerr << "IN: ReadMapper::write_mappings(string& file, bool filtered);" << endl;
		cerr << "\tERROR: could not open file= " << file << endl;
		return;
	}
	
	long int n_mappings;
	if (filtered) {
		n_mappings = _fmapping_ry.size();
		
		for (long int i = 0; i < n_mappings; ++i) {
			out << _fmapping_ry[i] << '\t' << _fmapping_ws[i] << '\t';
			out << _fmapping_mm[i] << '\t' << _fmapping_loc[i] << '\t';
			out << _fmapping_gene[i] << '\t' << _fmapping_gene_strand[i] << '\t';
			out << _fmapping_acession[i] << '\t' << _fmapping_org[i] << endl;
		}
		
	} else {
		n_mappings = _mapping_ry.size();
		
		for (long int i = 0; i < n_mappings; ++i) {
			out << _mapping_ry[i] << '\t' << _mapping_ws[i] << '\t';
			out << _mapping_mm[i] << '\t' << _mapping_loc[i] << '\t';
			out << _mapping_gene[i] << '\t' << _mapping_gene_strand[i] << '\t';
			out << _mapping_accession[i] << '\t' << _mapping_org[i] << endl;
		}
	}
}
//----------------------------- write_gene_hits(string& file) ----------------------------
//header: "org\taccession\tp=0,r=1\tgene_name\tgene_syn\tstrand\teHits\tcHits\tproduct\n"
void ReadMapper::write_gene_hits (string& file)
{
	ofstream outFH(file.c_str(), ios_base::app);
	if (!outFH) {
		cerr << "IN: ReadMapper::write_gene_hits(string &file);" << endl;
		cerr << "\tERROR: could not open file= " << file << endl;
	} else {
		// Go through the proteins
		for (int i=0; i < _n_proteins; ++i) {
			outFH << _seq_source << '\t' << _accession << '\t' << 0 << '\t';
			outFH << protein_names[i] << '\t' << protein_syns[i] << '\t';
			outFH << protein_strand[i] << '\t' << protein_eHits[i] << '\t';
			outFH << protein_cHits[i] << '\t' << protein_products[i] << endl;
		}
		
		//Go through the rnas
		for (int i = 0; i < _n_rna; ++i) {
			outFH << _seq_source << '\t' << _accession << '\t' << 1 << '\t';
			outFH << rna_names[i] << '\t' << rna_syns[i] << '\t';
			outFH << rna_strand[i] << '\t' << rna_eHits[i] << '\t';
			outFH << rna_cHits[i] << '\t' << rna_products[i] << endl;
		}
	}
}
//---------------------------- write_seq_hits(string& file) ------------------------------
// seq\taccession\ttotal\tclose\tfclose\texact\tfexact\teProt\tcProt\teRna\tcRna\t"
// eNCoding\tcNCoding\n
void ReadMapper::write_seq_hits (string& file)
{
	ofstream out(file.c_str(), ios_base::app);
	if (!out) {
		cerr << "IN: ReadMapper::write_seq_hits(string& file);" << endl;
		cerr << "\tERROR: could not open file= " << file << endl;
	} else {
		// Calc the total exact and close protein hits
		for (int i = 0; i < _n_proteins; ++i) {
			_n_exact_prot_hits += protein_eHits[i];
			_n_close_prot_hits += protein_cHits[i];
		}
		
		// Calc the total exact and close rna hits
		for (int i = 0; i < _n_rna; ++i) {
			_n_exact_rna_hits += rna_eHits[i];
			_n_close_rna_hits += rna_cHits[i];
		}
		
		// Calc the total exact and close non-coding hits
		long int n_fmappings = _fmapping_ry.size();
		for (long int i = 0; i < n_fmappings; ++i) {
			if (_fmapping_gene[i] == "-") {
				if (_fmapping_mm[i] == 0)
					++_n_eNCoding_hits;
				else
					++_n_cNCoding_hits;
			}
		}
		
		long int total = _n_fexact + _n_fclose;
		long int n_unfiltered_exact = _n_fexact + _n_exact_dup;
		
		out << _seq_source << '\t' << _accession << '\t';
		out << total << '\t' << _n_close << '\t' << _n_fclose << '\t';
		out << n_unfiltered_exact << '\t' << _n_fexact << '\t';
		out << _n_exact_prot_hits << '\t' << _n_close_prot_hits << '\t';
		out << _n_exact_rna_hits << '\t' << _n_close_rna_hits << '\t';
		out << _n_eNCoding_hits << '\t' << _n_cNCoding_hits << endl;
		
		out.close();
	}
}
